MOFA2 - downstream
Rasmus Rydbirk
2026-01-13
Please focus on association between factor 1 and striatum. I include calculations on data with and without imputation. Im glad to see they provide the same output, namely association between (subtype of) MSA and factor 1/striatum.
Setup
system("sudo apt-get update\nsudo apt-get install libnetcdf-dev -y")
library(magrittr)
library(dplyr)
library(MOFA2)
library(ggplot2)
library(UpSetR)
library(scHelper)
library(qs)
library(tidyr)
library(purrr)
library(tibble)
library(MOFAdata)
tt <- Sys.time()Helper functions
Weights with information about imputation
plot_top_weights_imp <- function (object, imp.list, view = 1, factors = 1, nfeatures = 10, abs = TRUE, scale = TRUE, sign = "all")
{
if (!is(object, "MOFA"))
stop("'object' has to be an instance of MOFA")
if (nfeatures <= 0)
stop("'nfeatures' has to be greater than 0")
if (sign == "all") {
abs <- TRUE
}
if (is.numeric(view))
view <- views_names(object)[view]
stopifnot(view %in% views_names(object))
view <- MOFA2:::.check_and_get_views(object, view)
factors <- MOFA2:::.check_and_get_factors(object, factors)
W <- get_weights(object, factors = factors, views = view,
as.data.frame = TRUE)
if (scale)
W$value <- W$value/max(abs(W$value))
W <- W[W$value != 0, ]
W$sign <- ifelse(W$value > 0, "+", "-")
if (sign == "positive") {
W <- W[W$value > 0, ]
}
else if (sign == "negative") {
W <- W[W$value < 0, ]
}
if (abs)
W$value <- abs(W$value)
W <- W[with(W, order(-abs(value))), ]
W <- as.data.frame(top_n(group_by(W, factor), n = nfeatures,
wt = value))
W$feature_id <- W$feature
if ((length(unique(W$view)) > 1) && (nfeatures > 0) && (any(duplicated(W[W$factor == factors[1], ]$feature_id)))) {
message("Duplicated feature names across views, we will add the view name as a prefix")
W$feature_id <- paste(W$view, W$feature, sep = "_")
}
W$feature_id <- factor(W$feature_id, levels = rev(unique(W$feature_id)))
# Add information on imputation
W %<>%
mutate(feature_stripped = as.character(feature_id) %>%
strsplit("_") %>%
sget(1))
W$imp <- imp.list[[view]] %>%
.[match(W$feature_stripped, names(.))] %>%
factor(levels = c(F, T), labels = c("Not imputed", "Imputed"))
se.assay <- SummarizedExperiment::assay(se.noimp) %>%
.[W$feature_stripped, grepl(view, colnames(.))]
assay.per.cond <- se.assay %>%
as.data.frame()
nsamples <- ncol(se.assay)
# Add per condition
assay_long <- se.assay %>%
as.data.frame() %>%
as_tibble(rownames = "protein") %>%
pivot_longer(-protein, names_to = "id", values_to = "value") %>%
left_join(meta %>% dplyr::select(id, Diagnosis), by = "id")
diag_counts <- assay_long %>%
group_by(protein, Diagnosis) %>%
summarise(
impn = sum(is.na(value)),
nsamp = n(),
impper = round(impn / nsamp * 100, 1),
.groups = "drop"
)
diag_wide <- diag_counts %>%
dplyr::select(protein, Diagnosis, impn, impper) %>%
pivot_wider(
id_cols = protein,
names_from = Diagnosis,
values_from = c(impn, impper),
names_sep = "_"
) %>%
.[match(rownames(se.assay), .$protein),] %>%
dplyr::select(-protein)
nsamples <- ncol(se.assay)
W %<>%
mutate(impn = se.assay %>%
apply(1, \(x) sum(is.na(x))) %>%
unname()) %>%
mutate(impper = round(impn / nsamples * 100, 1) %>%
unname()) %>%
cbind(diag_wide) %>%
dplyr::select(-feature, -view, -feature_stripped) %>%
mutate(feature_id = as.character(feature_id),
imp = as.character(imp),
factor = as.character(factor),
impper = as.character(impper)) %>%
{rbind(data.frame(factor = "Factor1", value = 0, sign = "-", feature_id = "Columns", imp = "Not imputed", impn = "Tot, n", impper = "Tot%", impn_CTRL = "CTRL, n", impn_MSA = "MSA, n", impn_PD = "PD, n", impn_PSP = "PSP, n", impper_CTRL = "CTRL%", impper_MSA = "MSA%", impper_PD = "PD%", impper_PSP = "PSP%"), .)} %>%
mutate(feature_id = factor(feature_id, levels = rev(unique(feature_id))))
if (any(is.na(imp.logi))) warning("There are NAs in the list of imputed proteins.")
# Plot
p <- ggplot(W, aes(x = feature_id, y = value, col = imp)) +
geom_point(size = 2) +
geom_segment(aes(xend = feature_id), linewidth = 0.75, yend = 0, show.legend = F) +
scale_discrete_manual(aesthetics = "col", values = c("Imputed" = "firebrick", "Not imputed" = "black")) +
coord_flip() +
labs(y = "Weight") +
theme_bw() +
theme(axis.title.x = element_text(color = "black"),
axis.title.y = element_blank(),
axis.text.y = element_text(size = rel(1.1), hjust = 1, color = "black"),
axis.text.x = element_text(color = "black"),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(),
legend.position = "right",
legend.title = element_blank(),
legend.text = element_text(color = "black"),
legend.key = element_rect(fill = "transparent"),
strip.text = element_text(size = rel(1.2)),
panel.background = element_blank(),
panel.spacing = unit(1, "lines"),
panel.grid.major.y = element_blank(),
) #+
# facet_wrap(~factor, nrow = 1, scales = "free")
if (sign == "negative")
p <- p + scale_x_discrete(position = "top")
if (abs) {
p <- p +
ylim(0, max(W$value) + 0.1) +
geom_text(label = W$sign,
y = max(W$value) + 0.1,
size = 10,
show.legend = F)
}
p <- p +
ylim(0, max(W$value) + 2) +
geom_text(label = W$impn,
y = max(W$value) + 0.2,
size = 3,
show.legend = F) +
geom_text(label = W$impper,
y = max(W$value) + 0.4,
size = 3,
show.legend = F) +
geom_text(label = W$impn_CTRL,
y = max(W$value) + 0.6,
size = 3,
show.legend = F) +
geom_text(label = W$impper_CTRL,
y = max(W$value) + 0.8,
size = 3,
show.legend = F) +
geom_text(label = W$impn_MSA,
y = max(W$value) + 1,
size = 3,
show.legend = F) +
geom_text(label = W$impper_MSA,
y = max(W$value) + 1.2,
size = 3,
show.legend = F) +
geom_text(label = W$impn_PD,
y = max(W$value) + 1.4,
size = 3,
show.legend = F) +
geom_text(label = W$impper_PD,
y = max(W$value) + 1.6,
size = 3,
show.legend = F) +
geom_text(label = W$impn_PSP,
y = max(W$value) + 1.8,
size = 3,
show.legend = F) +
geom_text(label = W$impper_PSP,
y = max(W$value) + 2,
size = 3,
show.legend = F)
return(p)
}Sample-wise contribution to views (areas)
Function
pcs <- function (object, samples = "all", group_by = NULL, return_data = FALSE, ...)
{
if (!is(object, "MOFA"))
stop("'object' has to be an instance of MOFA")
scores <- object@cache$contribution_scores %>%
reshape2::melt()
colnames(scores) <- c("sample", "view", "value")
if (is.null(group_by)) {
to.plot <- scores
if (return_data)
return(to.plot)
p <- ggplot(to.plot, aes(x = .data$view, y = .data$value)) +
geom_bar(aes(fill = view), stat = "identity", color = "black") +
facet_wrap(~sample) + labs(x = "", y = "Contribution score") +
theme_classic() + theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(), legend.position = "top",
legend.title = element_blank())
return(p)
}
else {
to.plot <- merge(scores, object@samples_metadata[, c("sample",
group_by)], by = "sample")
if (return_data)
return(to.plot)
p <- ggplot(to.plot, aes(x = .data$view, y = .data$value)) +
geom_boxplot(aes(fill = view)) + facet_wrap(as.formula(paste("~",
group_by))) + labs(x = "", y = "Contribution score") +
theme_classic() + theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(), legend.position = "top",
legend.title = element_blank())
return(p)
}
}Region-based, with NA, with imputation
https://raw.githack.com/bioFAM/MOFA2_tutorials/master/R_tutorials/downstream_analysis.html
Load model
MOFAobject.trained <- load_model("/work/proteomics/02_data/mofa_nogroups_withna_withimp.hdf5")## The value -2^31 was detected in the dataset.
## This has been converted to NA within R.## Warning in .quality_control(object, verbose = verbose): Factor(s) 1, 3 are strongly correlated with the total number of expressed features for at least one of your omics. Such factors appear when there are differences in the total 'levels' between your samples, *sometimes* because of poor normalisation in the preprocessing steps.Factors correlations
plot_factor_cor(MOFAobject.trained)
Explained variance per area
plot_variance_explained(MOFAobject.trained)
MOFAobject.trained@cache$variance_explained$r2_per_factor## $group1
## CER FC MED SN STR
## Factor1 2.339857817 2.41563320 2.753508091 2.576148510 24.133199453
## Factor2 0.067317486 5.63617349 7.398545742 7.600158453 0.005507469
## Factor3 0.015407801 0.23275614 0.007551908 13.743638992 0.852268934
## Factor4 11.814665794 0.04336834 0.039535761 0.004839897 2.210634947
## Factor5 4.032844305 2.28035450 0.046372414 2.125096321 5.322384834
## Factor6 0.014066696 2.31698751 8.772975206 1.247757673 0.294625759
## Factor7 1.019167900 0.17645359 1.351934671 5.049353838 1.239514351
## Factor8 2.159816027 0.01425147 2.627253532 1.513928175 2.204763889
## Factor9 0.418668985 1.20673180 0.087374449 0.857770443 4.263985157
## Factor10 0.006401539 0.75365901 0.007772446 0.007688999 4.487466812plot_variance_explained(MOFAobject.trained,
plot_total = T)[[2]]
calculate_variance_explained_per_sample(MOFAobject.trained)$group1 %>%
as.data.frame() %>%
mutate(diagnosis = MOFAobject.trained@samples_metadata$Diagnosis) %>%
reshape2::melt(id.vars = "diagnosis") %>%
ggplot(aes(variable, value, fill = diagnosis)) +
geom_boxplot() +
geom_point(position = position_jitterdodge(jitter.width = 0.05, dodge.width = 0.7)) +
theme_bw() +
theme(line = element_blank()) +
labs(x = "", y = "Variance explained", fill = "Diagnosis") +
scale_fill_manual(values = ggsci::pal_npg()(4))## Warning: Removed 27 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).
Covariate association
correlate_factors_with_covariates(MOFAobject.trained,
covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"),
plot="r",
return_data = F)## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...
correlate_factors_with_covariates(MOFAobject.trained,
covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"),
return_data = F,
alpha = 1)## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...
Factor values
plot_factor(MOFAobject.trained,
factor = 1:10,
color_by = "Subtype",
shape_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factor = 1:10,
color_by = "Subtype",
shape_by = "Diagnosis",
group_by = "Diagnosis",
dodge = T,
add_boxplot = T,
)
Factor 1
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "Diagnosis",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "Age",
shape_by = "Diagnosis",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "PMI",
shape_by = "Diagnosis",
dot_size = 3
)
Factor 2
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "Diagnosis",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "PMI",
shape_by = "Diagnosis",
dot_size = 3,
)
Factor 3
plot_factor(MOFAobject.trained,
factors = 3,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
Factor 4
plot_factor(MOFAobject.trained,
factors = 4,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
Factor 5
plot_factor(MOFAobject.trained,
factors = 5,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factors = 5,
color_by = "Diagnosis",
dot_size = 3
)
Factor 7
plot_factor(MOFAobject.trained,
factors = 7,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factors = 7,
color_by = "Diagnosis",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 7,
color_by = "Sex",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 7,
color_by = "Age",
dot_size = 3
)
Combined
plot_factors(MOFAobject.trained,
factors = 1:3,
color_by = "Diagnosis",
shape_by = "Subtype"
)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2## Warning: Groups with fewer than two data points have been dropped.## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf## Warning: Groups with fewer than two data points have been dropped.## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf## Warning: Groups with fewer than two data points have been dropped.## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
plot_factors(MOFAobject.trained,
factors = c(1,3),
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3) +
geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey") +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey")
plot_factors(MOFAobject.trained,
factors = c(1,2),
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3) +
geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey")
# geom_hline(yintercept = 1, linetype = "dashed", color = "grey")plot_factors(MOFAobject.trained,
factors = c(2,4),
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3) +
geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey") +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey")
Weights
se.imp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")$dat.clean
se.noimp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")$dat.norm
meta <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")$expDesign
areas <- names(se.imp@elementMetadata@listData)[grep("imputed", names(se.imp@elementMetadata@listData))] %>%
strsplit("_") %>%
sget(1)
imp.logi <- se.imp@elementMetadata@listData %>%
.[grepl("imputed", names(.))] %>%
setNames(areas)F1, STR
plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
view = "STR",
factor = 1,
nfeatures = 30,
scale = F
)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
plot_data_heatmap(MOFAobject.trained,
view = "STR", # view of interest
factor = 1, # factor of interest
features = 15, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = F,
annotation_samples = "Subtype",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "STARD10_STR",
group_by = "Subtype",
shape_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "UTF1_STR",
group_by = "Subtype",
shape_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "ELK4_STR",
group_by = "Subtype",
shape_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "STR",
factor = 1,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
F2
Shared proteins
comp <- c("FC", "MED", "SN") %>%
lapply(\(aa) {
tmp <- MOFAobject.trained@expectations$W[[aa]] %>%
as.data.frame() %>%
dplyr::select(Factor2) %>%
arrange(-abs(Factor2))
out <- list()
for (i in c(20, 50, 100)) {
out[[paste(aa, i, sep = "_")]] <- dplyr::slice(tmp, seq(i)) %>%
rownames() %>%
strsplit("_") %>%
sget(1)
}
return(out)
})
seq(3) %>%
lapply(\(x) lget(comp, x) %>%
setNames(c("FC", "MED", "SN")) %>%
fromList() %>%
upset(nsets = 3, title = paste0("Top ", c(20, 50, 100)[x])))## [[1]]
##
## [[2]]
##
## [[3]]
lget(comp, 3) %>%
setNames(c("FC", "MED", "SN")) %>%
unlist() %>%
table() %>%
data.frame() %>%
setNames(c("symbol", "Freq")) %>%
filter(Freq >= 2) %>%
arrange(-Freq, symbol)## symbol Freq
## 1 LIMS1 3
## 2 LRPPRC 3
## 3 SUCLA2 3
## 4 SUCLG1 3
## 5 ACOT8 2
## 6 ALG13 2
## 7 CHN1 2
## 8 COMMD2 2
## 9 CTSG 2
## 10 DCK 2
## 11 FSD1 2
## 12 GDPGP1 2
## 13 GLRX5 2
## 14 GLS 2
## 15 H2AX 2
## 16 H3-7 2
## 17 H4C4 2
## 18 HARS2 2
## 19 HLCS 2
## 20 IARS2 2
## 21 IMPACT 2
## 22 LARS2 2
## 23 LMCD1 2
## 24 MCCC1 2
## 25 NCK1 2
## 26 NCK2 2
## 27 NCKAP1L 2
## 28 NSUN2 2
## 29 NUDT4 2
## 30 NUDT9 2
## 31 PDE5A 2
## 32 PNPT1 2
## 33 POR 2
## 34 PTGR3 2
## 35 RELA 2
## 36 RPRD1B 2
## 37 SMAP2 2
## 38 STARD10 2
## 39 WDSUB1 2FC
plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
view = "FC",
factor = 2,
nfeatures = 30,
scale = F
)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
plot_data_heatmap(MOFAobject.trained,
view = "FC", # view of interest
factor = 2, # factor of interest
features = 20, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "LRPPRC_FC",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "SUCLA2_FC",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "PRTN3_FC",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "FC",
factor = 2,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
MED
plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
view = "MED",
factor = 2,
nfeatures = 30,
scale = F
)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
plot_data_heatmap(MOFAobject.trained,
view = "MED", # view of interest
factor = 2, # factor of interest
features = 20, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "TMEM119_MED",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "SUCLA2_MED",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "MED",
factor = 2,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
SN
plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
view = "SN",
factor = 2,
nfeatures = 82,
scale = F
)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
plot_data_heatmap(MOFAobject.trained,
view = "SN", # view of interest
factor = 2, # factor of interest
features = 30, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "STARD10_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "RPL34_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "SN",
factor = 2,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
F3, SN
plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
view = "SN",
factor = 3,
nfeatures = 30,
scale = F
)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
plot_data_heatmap(MOFAobject.trained,
view = "SN", # view of interest
factor = 3, # factor of interest
features = 20, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 3,
color_by = "RPL32_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 3,
color_by = "PDE1B_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "SN",
factor = 3,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
tSNE
model <- run_tsne(MOFAobject.trained, perplexity = 10)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3
)## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the MOFA2 package.
## Please report the issue at <https://github.com/bioFAM/MOFA2>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Diagnosis",
# shape_by = "Diagnosis",
dot_size = 3
)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Factor1",
shape_by = "Diagnosis",
dot_size = 3
)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Factor2",
shape_by = "Diagnosis",
dot_size = 3
)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Factor3",
shape_by = "Diagnosis",
dot_size = 3
)
Sample-wise contribution scores
MOFAobject.trained <- calculate_contribution_scores(object = MOFAobject.trained)
samples_metadata(MOFAobject.trained)## Age Diagnosis Duration Hemisphere Origin PMI Sex Subtype group sample
## 1 65 PSP 7 Left BBH 28.0 Male RS group1 1370
## 2 73 MSA 7 Left BBH 24.0 Male P group1 1379
## 3 74 MSA 4 Left BBH 56.0 Female C group1 1388
## 4 59 MSA 8 Right BBH 24.0 Female C group1 1391
## 5 56 MSA 8 Right BBH 72.0 Male C+P group1 1392
## 6 80 PSP 9 Right BBH 36.0 Male RS group1 1394
## 7 55 MSA 3 Left BBH 31.0 Female C group1 1398
## 8 66 PD 13 Right BBH 53.0 Male None group1 1400
## 9 69 MSA 8 Left BBH 24.0 Female P group1 1406
## 10 56 MSA 7 Right BBH 58.0 Male P group1 1407
## 11 70 MSA 9 Right BBH 26.0 Male P group1 1436
## 12 79 PD 10 Left BBH 118.0 Male None group1 1446
## 13 68 MSA 7 Right BBH 70.0 Male P group1 1462
## 14 91 CTRL NA Left BBH 24.0 Female None group1 1464
## 15 68 CTRL NA Right BBH 24.0 Male None group1 1467
## 16 91 CTRL NA Left BBH 24.0 Female None group1 1490
## 17 82 CTRL NA Right BBH 48.0 Male None group1 1492
## 18 86 CTRL NA Left BBH 45.0 Female None group1 1494
## 19 97 CTRL NA Right BBH 90.0 Female None group1 1495
## 20 82 PSP 5 Left BBH 24.0 Male None group1 1385
## 21 87 PSP 8 Right BBH 30.0 Female None group1 1387
## 22 67 PSP 9 Left BBH 30.0 Male None group1 1397
## 23 71 PSP 4 Right BBH 30.0 Male RS group1 1401
## 24 69 PSP 9 Right BBH 29.0 Male Cortical group1 1410
## 25 81 PD 9 Left BBH 34.0 Male None group1 1411
## 26 70 PSP 6 Left BBH 19.0 Male RS group1 1414
## 27 71 PSP 8 Right BBH 101.0 Male Cortical group1 1420
## 28 64 PD 13 Left BBH 117.5 Male None group1 1440
## 29 88 CTRL NA Right BBH 24.0 Female None group1 1453
## 30 72 PD 6 Right BBH 89.0 Male None group1 1457
## 31 84 CTRL NA Right BBH 96.0 Female None group1 1466
## 32 75 PSP 7 Right BBH 15.0 Female RS group1 1381
## 33 95 PD 7 Left BBH 86.0 Female None group1 1412
## 34 91 PD 10 Left BBH 39.0 Female None group1 1458
## 35 86 PD 12 Left BBH 39.0 Male None group1 1424
## 36 75 PD NA Left BBH 48.0 Male None group1 1449
## 37 48 CTRL NA Whole BBH 48.0 Female None group1 1461
## 38 91 CTRL NA Left BBH 96.0 Female None group1 1465
## 39 75 PD NA Left BBH 32.0 Female None group1 1429
## 40 62 PD NA Right BBH 34.0 Male None group1 1469
## CER_contribution FC_contribution MED_contribution SN_contribution
## 1 0.15082781 0.06216384 0.12712092 0.46087914
## 2 0.14714621 0.27634369 0.17855782 0.17427280
## 3 0.12335640 0.22708067 0.21952391 0.29500351
## 4 NaN NaN NaN NaN
## 5 NaN NaN NaN NaN
## 6 0.22138774 0.13635375 0.23508441 0.24971132
## 7 NaN NaN NaN NaN
## 8 0.23427308 0.24046448 0.20627061 0.16327756
## 9 0.23135490 0.21387768 0.14590016 0.16719717
## 10 NaN NaN NaN NaN
## 11 NaN NaN NaN NaN
## 12 0.21281771 0.23732481 0.14386542 0.22032829
## 13 NaN NaN NaN NaN
## 14 0.26906764 0.25865746 0.16421821 0.13261004
## 15 NaN NaN NaN NaN
## 16 NaN NaN NaN NaN
## 17 0.13177486 0.16031751 0.03116276 0.49728780
## 18 0.22941892 0.22196431 0.05917558 0.31262607
## 19 0.18399271 0.20715479 0.27007271 0.19715939
## 20 NaN NaN NaN NaN
## 21 0.25899102 0.00000000 0.14081864 0.16015999
## 22 0.32527763 0.11095035 0.08469555 0.27869658
## 23 0.26404507 0.10151794 0.10497410 0.33232702
## 24 NaN NaN NaN NaN
## 25 0.09431939 0.01175833 0.10577559 0.63934940
## 26 0.27749772 0.07444969 0.26094747 0.13652442
## 27 0.37322884 0.00000000 0.23567577 0.14646064
## 28 0.10797343 0.32450990 0.27192454 0.23417692
## 29 0.16059805 0.08809670 0.28351085 0.24027018
## 30 0.21043819 0.32227562 0.15072773 0.15046156
## 31 0.35634937 0.15094665 0.22551296 0.09514165
## 32 NaN NaN NaN NaN
## 33 NaN NaN NaN NaN
## 34 NaN NaN NaN NaN
## 35 NaN NaN NaN NaN
## 36 NaN NaN NaN NaN
## 37 NaN NaN NaN NaN
## 38 NaN NaN NaN NaN
## 39 NaN NaN NaN NaN
## 40 NaN NaN NaN NaN
## STR_contribution
## 1 0.19900829
## 2 0.22367948
## 3 0.13503551
## 4 NaN
## 5 NaN
## 6 0.15746278
## 7 NaN
## 8 0.15571427
## 9 0.24167008
## 10 NaN
## 11 NaN
## 12 0.18566378
## 13 NaN
## 14 0.17544665
## 15 NaN
## 16 NaN
## 17 0.17945707
## 18 0.17681512
## 19 0.14162040
## 20 NaN
## 21 0.44003035
## 22 0.20037988
## 23 0.19713586
## 24 NaN
## 25 0.14879729
## 26 0.25058070
## 27 0.24463476
## 28 0.06141521
## 29 0.22752422
## 30 0.16609690
## 31 0.17204937
## 32 NaN
## 33 NaN
## 34 NaN
## 35 NaN
## 36 NaN
## 37 NaN
## 38 NaN
## 39 NaN
## 40 NaNpcs(MOFAobject.trained, group_by = "group")## Warning: Removed 95 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
pcs(MOFAobject.trained)## Warning: Removed 95 rows containing missing values or values outside the scale range
## (`geom_bar()`).
Region-based, with NA, no imputation
https://raw.githack.com/bioFAM/MOFA2_tutorials/master/R_tutorials/downstream_analysis.html
Load model
MOFAobject.trained <- load_model("/work/proteomics/02_data/mofa_nogroups_withna_noimp.hdf5")## The value -2^31 was detected in the dataset.
## This has been converted to NA within R.## Warning in .quality_control(object, verbose = verbose): Factor(s) 1, 3 are strongly correlated with the total number of expressed features for at least one of your omics. Such factors appear when there are differences in the total 'levels' between your samples, *sometimes* because of poor normalisation in the preprocessing steps.Factors correlations
plot_factor_cor(MOFAobject.trained)
# =========================
# Setup
# =========================
W_list <- MOFAobject.trained@expectations$W # list of view matrices (features x factors)
# Fixed view order (change if you prefer a different order)
view_order <- c("CER", "FC", "MED", "SN", "STR")
W_list <- W_list[view_order]
# Ensure factor names exist and are consistent (Factor1..FactorN)
normalize_factor_names <- function(m) {
if (is.null(colnames(m)) || any(nchar(colnames(m)) == 0)) {
colnames(m) <- paste0("Factor", seq_len(ncol(m)))
}
m
}
W_list <- lapply(W_list, normalize_factor_names)
# Helper: strip trailing "_VIEW" suffix (e.g., IGKV3-7_CER -> IGKV3-7)
strip_suffix <- function(ids) sub("_[^_]+$", "", ids)
# Pairwise correlations of W (features x factors) between two views on common features
pair_cor_views <- function(W_i, W_j, method = "pearson") {
rid_i <- strip_suffix(rownames(W_i))
rid_j <- strip_suffix(rownames(W_j))
common <- intersect(rid_i, rid_j)
# If overlap is tiny, return NA matrix with proper dimnames
if (length(common) < 3) {
out <- matrix(NA_real_, nrow = ncol(W_i), ncol = ncol(W_j))
rownames(out) <- colnames(W_i)
colnames(out) <- colnames(W_j)
return(out)
}
A <- W_i[match(common, rid_i), , drop = FALSE]
B <- W_j[match(common, rid_j), , drop = FALSE]
out <- cor(A, B, use = "pairwise.complete.obs", method = method)
# Preserve factor names
rownames(out) <- colnames(W_i)
colnames(out) <- colnames(W_j)
out
}
# =========================
# Compute only upper triangle (including diagonal) panels
# =========================
method <- "pearson" # or "spearman"
top_k <- 3 # number of top correlations to annotate per panel
pairwise_mats <- list()
for (i in seq_along(view_order)) {
for (j in seq_along(view_order)) {
if (i <= j) { # upper triangle + diagonal
v1 <- view_order[i]; v2 <- view_order[j]
lbl <- paste(v1, v2, sep = " vs ")
pairwise_mats[[lbl]] <- pair_cor_views(W_list[[v1]], W_list[[v2]], method = method)
}
}
}
# =========================
# Melt for ggplot facets
# =========================
if (!requireNamespace("reshape2", quietly = TRUE)) install.packages("reshape2")
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
melt_pair <- function(mat, v1, v2) {
df <- reshape2::melt(mat)
colnames(df) <- c("Factor1", "Factor2", "cor")
df$view1 <- v1
df$view2 <- v2
df
}
df_all <- do.call(
rbind,
lapply(names(pairwise_mats), function(lbl) {
parts <- strsplit(lbl, " vs ", fixed = TRUE)[[1]]
melt_pair(pairwise_mats[[lbl]], parts[1], parts[2])
})
)
# Factor ordering (ensures axes read Factor1..FactorN everywhere)
all_factors <- sort(unique(unlist(lapply(W_list, colnames))))
df_all$Factor1 <- factor(df_all$Factor1, levels = all_factors)
df_all$Factor2 <- factor(df_all$Factor2, levels = all_factors)
# =========================
# Build top-k overlay (per panel) EXCLUDING r == 1 and diagonal self-pairs
# =========================
split_panels <- split(df_all, interaction(df_all$view1, df_all$view2, drop = TRUE))
top_list <- lapply(split_panels, function(d) {
# Drop NA and exact 1.0 correlations
d <- subset(d, !is.na(cor) & cor != 1)
# In within-view panels, exclude diagonal cells (Factor1 == Factor2)
# which would be 1 by definition; already excluded by cor != 1,
# but keep the explicit filter in case of numerical jitter.
if (nrow(d) > 0 && unique(d$view1) == unique(d$view2)) {
d <- subset(d, as.character(Factor1) != as.character(Factor2))
}
# Order by absolute correlation descending
d <- d[order(-abs(d$cor)), , drop = FALSE]
# Keep top-k
if (nrow(d) == 0) return(d)
d[seq_len(min(nrow(d), top_k)), , drop = FALSE]
})
top_df <- do.call(rbind, top_list)
# =========================
# Plot
# =========================
labeller_fn <- ggplot2::labeller(
view1 = function(x) paste0("View 1: ", x),
view2 = function(x) paste0("View 2: ", x)
)
p <- ggplot2::ggplot(df_all, ggplot2::aes(x = Factor2, y = Factor1, fill = cor)) +
ggplot2::geom_tile() +
# Overlay: mark top-k tiles with a ring + correlation value
ggplot2::geom_point(
data = top_df,
shape = 21, stroke = 0.9, size = 3.8, colour = "black", fill = NA
) +
ggplot2::geom_text(
data = top_df,
ggplot2::aes(label = sprintf("%.2f", cor)),
size = 3, fontface = "bold", colour = "black"
) +
ggplot2::scale_fill_gradient2(
low = "navy", mid = "white", high = "firebrick3",
midpoint = 0, limits = c(-1, 1), oob = scales::squish, name = "r"
) +
ggplot2::facet_grid(view1 ~ view2, labeller = labeller_fn, drop = TRUE) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(
strip.text = ggplot2::element_text(face = "bold"),
axis.text.x = ggplot2::element_text(angle = 45, hjust = 1),
panel.spacing = grid::unit(0.6, "lines"),
legend.position = "right"
) +
ggplot2::labs(
title = sprintf("MOFA W: Factor × Factor correlations (upper triangle), method = %s", method),
subtitle = sprintf("Top-%d factor pairs per panel highlighted (excluding r = 1); panels use pairwise intersections of features", top_k),
x = "Factors in View 2",
y = "Factors in View 1"
)
print(p)
Explained variance per area
plot_variance_explained(MOFAobject.trained)
MOFAobject.trained@cache$variance_explained$r2_per_factor## $group1
## CER FC MED SN STR
## Factor1 0.072622299 3.231221437 2.674627304 0.01284480 36.455994844
## Factor2 14.094030857 3.056883812 9.029161930 1.87025070 5.805820227
## Factor3 6.209421158 3.535330296 0.028491020 19.61889863 1.466298103
## Factor4 6.710535288 5.907905102 0.435268879 2.19438076 10.662490129
## Factor5 0.005841255 1.038879156 17.748779058 0.00308156 4.222744703
## Factor6 2.574747801 1.141369343 6.447696686 11.38798594 0.010675192
## Factor7 7.396245003 6.605488062 0.034427643 7.38848448 0.007265806
## Factor8 0.006121397 0.005877018 0.008529425 10.15772820 0.357550383plot_variance_explained(MOFAobject.trained,
plot_total = T)[[2]]
Covariate association
correlate_factors_with_covariates(MOFAobject.trained,
covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"),
plot="r",
return_data = F)## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...
correlate_factors_with_covariates(MOFAobject.trained,
covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"),
return_data = F,
alpha = 1)## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...
Factor values
plot_factor(MOFAobject.trained,
factor = 1:8,
color_by = "Subtype",
shape_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factor = 1:8,
color_by = "Subtype",
shape_by = "Diagnosis",
group_by = "Diagnosis",
dodge = T,
add_boxplot = T,
)
Factor 1
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "Diagnosis",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "Age",
shape_by = "Diagnosis",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "PMI",
shape_by = "Diagnosis",
dot_size = 3
)
Factor 2
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "Diagnosis",
dot_size = 3
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "Sex",
shape_by = "Diagnosis",
dot_size = 3,
)
Factor 3
plot_factor(MOFAobject.trained,
factors = 3,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
Factor 4
plot_factor(MOFAobject.trained,
factors = 4,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
Factor 5
plot_factor(MOFAobject.trained,
factors = 5,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis"
)
Factor 6
plot_factor(MOFAobject.trained,
factors = 6,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis",
)
plot_factor(MOFAobject.trained,
factors = 6,
color_by = "Diagnosis",
dot_size = 3
)
Factor 8
plot_factor(MOFAobject.trained,
factors = 8,
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3,
group_by = "Diagnosis",
)
plot_factor(MOFAobject.trained,
factors = 8,
color_by = "Diagnosis",
dot_size = 3
)
Combined
plot_factors(MOFAobject.trained,
factors = c(1,3),
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3) +
geom_vline(xintercept = 0.9, linetype = "dashed", color = "grey") +
geom_hline(yintercept = 0.9, linetype = "dashed", color = "grey")
plot_factors(MOFAobject.trained,
factors = c(1,2),
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3) +
geom_vline(xintercept = 1, linetype = "dashed", color = "grey") +
geom_hline(yintercept = -1.1, linetype = "dashed", color = "grey")
plot_factors(MOFAobject.trained,
factors = c(6,8),
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3) +
geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey")
# geom_hline(yintercept = 1, linetype = "dashed", color = "grey")Weights
F1, STR
plot_top_weights(MOFAobject.trained,
view = "STR",
factor = 1,
nfeatures = 30,
scale = F
)
plot_data_heatmap(MOFAobject.trained,
view = "STR", # view of interest
factor = 1, # factor of interest
features = 15, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = F,
annotation_samples = "Subtype",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "PDE1B_STR",
group_by = "Subtype",
shape_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "UTF1_STR",
group_by = "Subtype",
shape_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 1,
color_by = "ATCAY_STR",
group_by = "Subtype",
shape_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "STR",
factor = 1,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
F2
Shared proteins
comp <- c("FC", "MED", "SN") %>%
lapply(\(aa) {
tmp <- MOFAobject.trained@expectations$W[[aa]] %>%
as.data.frame() %>%
dplyr::select(Factor2) %>%
arrange(-abs(Factor2))
out <- list()
for (i in c(20, 50, 100)) {
out[[paste(aa, i, sep = "_")]] <- dplyr::slice(tmp, seq(i)) %>%
rownames() %>%
strsplit("_") %>%
sget(1)
}
return(out)
})
seq(3) %>%
lapply(\(x) lget(comp, x) %>%
setNames(c("FC", "MED", "SN")) %>%
fromList() %>%
upset(nsets = 3, title = paste0("Top ", c(20, 50, 100)[x])))## [[1]]
##
## [[2]]
##
## [[3]]
FC
plot_top_weights(MOFAobject.trained,
view = "FC",
factor = 2,
nfeatures = 30,
scale = F
)
plot_data_heatmap(MOFAobject.trained,
view = "FC", # view of interest
factor = 2, # factor of interest
features = 20, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "PTGR3_FC",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "OGDHL_FC",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "GRPEL1_FC",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "FC",
factor = 2,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
MED
plot_top_weights(MOFAobject.trained,
view = "MED",
factor = 2,
nfeatures = 30,
scale = F
)
plot_data_heatmap(MOFAobject.trained,
view = "MED", # view of interest
factor = 2, # factor of interest
features = 20, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "FSD1_MED",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "SUCLA2_MED",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "MED",
factor = 2,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
SN
plot_top_weights(MOFAobject.trained,
view = "SN",
factor = 2,
nfeatures = 30,
scale = F
)
plot_data_heatmap(MOFAobject.trained,
view = "SN", # view of interest
factor = 2, # factor of interest
features = 30, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "CBR1_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 2,
color_by = "TRAP1_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "SN",
factor = 2,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
F3, SN
plot_top_weights(MOFAobject.trained,
view = "SN",
factor = 3,
nfeatures = 30,
scale = F
)
plot_data_heatmap(MOFAobject.trained,
view = "SN", # view of interest
factor = 3, # factor of interest
features = 20, # number of features to plot (they are selected by weight)
cluster_rows = T,
cluster_cols = T,
show_rownames = TRUE,
show_colnames = T,
annotation_samples = "Diagnosis",
denoise = T
)
plot_factor(MOFAobject.trained,
factors = 3,
color_by = "XPO1_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_factor(MOFAobject.trained,
factors = 3,
color_by = "CSE1L_SN",
group_by = "Diagnosis",
add_violin = TRUE,
)
plot_data_scatter(MOFAobject.trained,
view = "SN",
factor = 3,
features = 4,
color_by = "Subtype",
shape_by = "Diagnosis"
)
tSNE
model <- run_tsne(MOFAobject.trained, perplexity = 10)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Subtype",
shape_by = "Diagnosis",
dot_size = 3
)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Diagnosis",
# shape_by = "Diagnosis",
dot_size = 3
)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Factor1",
shape_by = "Diagnosis",
dot_size = 3
)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Factor2",
shape_by = "Diagnosis",
dot_size = 3
)
plot_dimred(model,
method = "TSNE", # method can be either "TSNE" or "UMAP"
color_by = "Factor3",
shape_by = "Diagnosis",
dot_size = 3
)
Sample-wise contribution scores
MOFAobject.trained <- calculate_contribution_scores(object = MOFAobject.trained)
samples_metadata(MOFAobject.trained)## Age Diagnosis Duration Hemisphere Origin PMI Sex Subtype group sample
## 1 65 PSP 7 Left BBH 28.0 Male RS group1 1370
## 2 73 MSA 7 Left BBH 24.0 Male P group1 1379
## 3 74 MSA 4 Left BBH 56.0 Female C group1 1388
## 4 59 MSA 8 Right BBH 24.0 Female C group1 1391
## 5 56 MSA 8 Right BBH 72.0 Male C+P group1 1392
## 6 80 PSP 9 Right BBH 36.0 Male RS group1 1394
## 7 55 MSA 3 Left BBH 31.0 Female C group1 1398
## 8 66 PD 13 Right BBH 53.0 Male None group1 1400
## 9 69 MSA 8 Left BBH 24.0 Female P group1 1406
## 10 56 MSA 7 Right BBH 58.0 Male P group1 1407
## 11 70 MSA 9 Right BBH 26.0 Male P group1 1436
## 12 79 PD 10 Left BBH 118.0 Male None group1 1446
## 13 68 MSA 7 Right BBH 70.0 Male P group1 1462
## 14 91 CTRL NA Left BBH 24.0 Female None group1 1464
## 15 68 CTRL NA Right BBH 24.0 Male None group1 1467
## 16 91 CTRL NA Left BBH 24.0 Female None group1 1490
## 17 82 CTRL NA Right BBH 48.0 Male None group1 1492
## 18 86 CTRL NA Left BBH 45.0 Female None group1 1494
## 19 97 CTRL NA Right BBH 90.0 Female None group1 1495
## 20 82 PSP 5 Left BBH 24.0 Male None group1 1385
## 21 87 PSP 8 Right BBH 30.0 Female None group1 1387
## 22 67 PSP 9 Left BBH 30.0 Male None group1 1397
## 23 71 PSP 4 Right BBH 30.0 Male RS group1 1401
## 24 69 PSP 9 Right BBH 29.0 Male Cortical group1 1410
## 25 81 PD 9 Left BBH 34.0 Male None group1 1411
## 26 70 PSP 6 Left BBH 19.0 Male RS group1 1414
## 27 71 PSP 8 Right BBH 101.0 Male Cortical group1 1420
## 28 64 PD 13 Left BBH 117.5 Male None group1 1440
## 29 88 CTRL NA Right BBH 24.0 Female None group1 1453
## 30 72 PD 6 Right BBH 89.0 Male None group1 1457
## 31 84 CTRL NA Right BBH 96.0 Female None group1 1466
## 32 75 PSP 7 Right BBH 15.0 Female RS group1 1381
## 33 95 PD 7 Left BBH 86.0 Female None group1 1412
## 34 91 PD 10 Left BBH 39.0 Female None group1 1458
## 35 86 PD 12 Left BBH 39.0 Male None group1 1424
## 36 75 PD NA Left BBH 48.0 Male None group1 1449
## 37 48 CTRL NA Whole BBH 48.0 Female None group1 1461
## 38 91 CTRL NA Left BBH 96.0 Female None group1 1465
## 39 75 PD NA Left BBH 32.0 Female None group1 1429
## 40 62 PD NA Right BBH 34.0 Male None group1 1469
## CER_contribution FC_contribution MED_contribution SN_contribution
## 1 0.14214162 0.07123318 0.2335368 0.3917811
## 2 0.06725632 0.31186582 0.2231970 0.1263306
## 3 0.18876339 0.21377168 0.2352642 0.2225193
## 4 NaN NaN NaN NaN
## 5 NaN NaN NaN NaN
## 6 0.16822988 0.23762093 0.1425290 0.2873132
## 7 NaN NaN NaN NaN
## 8 0.22864984 0.23697852 0.2687660 0.1936276
## 9 0.14957799 0.11660456 0.2413496 0.1457976
## 10 NaN NaN NaN NaN
## 11 NaN NaN NaN NaN
## 12 0.21213553 0.26038667 0.1666308 0.1829730
## 13 NaN NaN NaN NaN
## 14 0.26612438 0.20493846 0.2456342 0.1432734
## 15 NaN NaN NaN NaN
## 16 NaN NaN NaN NaN
## 17 0.17166178 0.14881793 0.2031213 0.3280382
## 18 0.16583917 0.16992017 0.0645116 0.4507353
## 19 0.18931714 0.27299803 0.1633389 0.1954027
## 20 NaN NaN NaN NaN
## 21 0.25027903 0.07958009 0.1626316 0.3382785
## 22 0.29430106 0.04239646 0.1040508 0.3707292
## 23 0.13338611 0.16168388 0.1905265 0.2708038
## 24 NaN NaN NaN NaN
## 25 0.18935370 0.00000000 0.1647600 0.5344074
## 26 0.27797464 0.00000000 0.2874644 0.2305561
## 27 0.35424787 0.02548541 0.1898440 0.2271896
## 28 0.19559509 0.24869516 0.2153742 0.2067256
## 29 0.15693131 0.05634121 0.3298812 0.2506410
## 30 0.21998408 0.28382531 0.1515171 0.1386265
## 31 0.14490946 0.20805933 0.3463598 0.1873873
## 32 NaN NaN NaN NaN
## 33 NaN NaN NaN NaN
## 34 NaN NaN NaN NaN
## 35 NaN NaN NaN NaN
## 36 NaN NaN NaN NaN
## 37 NaN NaN NaN NaN
## 38 NaN NaN NaN NaN
## 39 NaN NaN NaN NaN
## 40 NaN NaN NaN NaN
## STR_contribution
## 1 0.16130730
## 2 0.27135021
## 3 0.13968144
## 4 NaN
## 5 NaN
## 6 0.16430709
## 7 NaN
## 8 0.07197804
## 9 0.34667028
## 10 NaN
## 11 NaN
## 12 0.17787409
## 13 NaN
## 14 0.14002964
## 15 NaN
## 16 NaN
## 17 0.14836080
## 18 0.14899375
## 19 0.17894316
## 20 NaN
## 21 0.16923073
## 22 0.18852240
## 23 0.24359972
## 24 NaN
## 25 0.11147885
## 26 0.20400486
## 27 0.20323313
## 28 0.13360999
## 29 0.20620532
## 30 0.20604697
## 31 0.11328414
## 32 NaN
## 33 NaN
## 34 NaN
## 35 NaN
## 36 NaN
## 37 NaN
## 38 NaN
## 39 NaN
## 40 NaNpcs(MOFAobject.trained, group_by = "group")## Warning: Removed 95 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
pcs(MOFAobject.trained)## Warning: Removed 95 rows containing missing values or values outside the scale range
## (`geom_bar()`).
Session info
tt - Sys.time()## Time difference of -1.257868 minssessionInfo()## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
##
## Matrix products: default
## BLAS/LAPACK: /opt/intel/oneapi/mkl/2025.3/lib/libmkl_gf_lp64.so.2; LAPACK version 3.12.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Berlin
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MOFAdata_1.24.0 tibble_3.2.1 purrr_1.0.4 tidyr_1.3.1
## [5] qs_0.27.3 scHelper_0.0.5 UpSetR_1.4.0 ggplot2_3.5.2
## [9] MOFA2_1.18.0 dplyr_1.1.4 magrittr_2.0.3
##
## loaded via a namespace (and not attached):
## [1] mnormt_2.1.1 gridExtra_2.3
## [3] rlang_1.1.6 clue_0.3-66
## [5] GetoptLong_1.0.5 matrixStats_1.5.0
## [7] compiler_4.5.2 mgcv_1.9-3
## [9] sccore_1.0.5 dir.expiry_1.16.0
## [11] png_0.1-8 vctrs_0.6.5
## [13] reshape2_1.4.4 stringr_1.5.1
## [15] pkgconfig_2.0.3 shape_1.4.6.1
## [17] crayon_1.5.3 fastmap_1.2.0
## [19] backports_1.5.0 XVector_0.48.0
## [21] labeling_0.4.3 rmarkdown_2.29
## [23] UCSC.utils_1.4.0 xfun_0.52
## [25] cachem_1.1.0 GenomeInfoDb_1.44.0
## [27] jsonlite_2.0.0 rhdf5filters_1.20.0
## [29] DelayedArray_0.34.1 Rhdf5lib_1.30.0
## [31] psych_2.5.3 broom_1.0.8
## [33] parallel_4.5.2 cluster_2.1.8.1
## [35] R6_2.6.1 bslib_0.9.0
## [37] stringi_1.8.7 RColorBrewer_1.1-3
## [39] reticulate_1.42.0 GGally_2.2.1
## [41] car_3.1-3 leidenAlg_1.1.5
## [43] GenomicRanges_1.60.0 jquerylib_0.1.4
## [45] Rcpp_1.0.14 SummarizedExperiment_1.38.1
## [47] iterators_1.0.14 knitr_1.50
## [49] IRanges_2.42.0 splines_4.5.2
## [51] Matrix_1.7-3 igraph_2.1.4
## [53] tidyselect_1.2.1 rstudioapi_0.17.1
## [55] dichromat_2.0-0.1 abind_1.4-8
## [57] yaml_2.3.10 stringfish_0.16.0
## [59] doParallel_1.0.17 codetools_0.2-20
## [61] lattice_0.22-7 plyr_1.8.9
## [63] Biobase_2.68.0 basilisk.utils_1.20.0
## [65] withr_3.0.2 evaluate_1.0.3
## [67] Rtsne_0.17 ggstats_0.9.0
## [69] RcppParallel_5.1.10 circlize_0.4.16
## [71] ggpubr_0.6.0 pillar_1.10.2
## [73] filelock_1.0.3 carData_3.0-5
## [75] MatrixGenerics_1.20.0 corrplot_0.95
## [77] foreach_1.5.2 stats4_4.5.2
## [79] generics_0.1.4 S4Vectors_0.46.0
## [81] scales_1.4.0 RApiSerialize_0.1.4
## [83] glue_1.8.0 pheatmap_1.0.12
## [85] tools_4.5.2 ggsignif_0.6.4
## [87] forcats_1.0.1 cowplot_1.1.3
## [89] rhdf5_2.52.0 grid_4.5.2
## [91] conos_1.5.2 colorspace_2.1-1
## [93] GenomeInfoDbData_1.2.14 nlme_3.1-168
## [95] basilisk_1.20.0 HDF5Array_1.36.0
## [97] Formula_1.2-5 cli_3.6.5
## [99] S4Arrays_1.8.0 ComplexHeatmap_2.24.0
## [101] uwot_0.2.3 gtable_0.3.6
## [103] rstatix_0.7.2 ggsci_3.2.0
## [105] sass_0.4.10 digest_0.6.37
## [107] BiocGenerics_0.54.0 SparseArray_1.8.0
## [109] ggrepel_0.9.6 rjson_0.2.23
## [111] farver_2.1.2 htmltools_0.5.8.1
## [113] lifecycle_1.0.4 httr_1.4.7
## [115] h5mread_1.0.0 GlobalOptions_0.1.2